/* 2010: Copyright by Massimo Bernaschi, Giorgio Parisi, Leonardo Parisi */
#include "cudamacro.h"
#include "mpiexch.h"

#include "hsgAMGPUGM.h"
#include "timing.h"

#include "dictionary.h"
#include "iniparser.h"
#include "util.h"
#if defined(DEBUG)
#include "util/cuPrintf.cu"
#endif

#include "MersenneTwister.h"
#if !defined(HOSTREDU)
#include "threadFenceReduction_kernel.cu"
#endif

extern void StartMpi(int *, int *, int *, char ***);
extern struct mpiexch *SetUpMpi(int, int, int);
extern void ExchMpi(struct mpiexch*);
extern REALENE SumEnergy(REALENE, int);
extern void StopMpi(void);
//-------------------------------------------------------------------------------------------------------
// Usage
//-------------------------------------------------------------------------------------------------------
void Usage(char *program) {
    fprintf(stderr,"Usage: %s -i inputfile [-r spin_file] [-o (over-relax only)] [-b (heat bath only)]\n",program);
}


//-------------------------------------------------------------------------------------------------------
// Init
//-------------------------------------------------------------------------------------------------------
void Init(void){
	
    if(sizeof(REAL)==sizeof(float)) {
        zero=0.0f;
        one=1.0f;
        two=2.0f;
        three=3.0f;
        zerofive=0.5f;
        zerotwentyfive=0.25f;
    } else if(sizeof(REAL)==sizeof(double)) {
        zero=0.0;
        one=1.0;
        two=2.0;
        three=3.0;
        zerofive=0.5;
        zerotwentyfive=0.25;
    } else {
        writelog(TRUE,APPLICATION_RC,"Invalid size of REAL: %d",sizeof(REAL));
    }
	
    zeroEne = 0.0;
	
    h_Ra = (REAL*) Malloc(ms_s);
    h_Rb = (REAL*) Malloc(ms_s);
    h_Rc = (REAL*) Malloc(ms_s);
    h_Ba = (REAL*) Malloc(ms_s);
    h_Bb = (REAL*) Malloc(ms_s);
    h_Bc = (REAL*) Malloc(ms_s);
	
    h_Jpx = (REAL*) Malloc(ms_j);
    h_Jpy = (REAL*) Malloc(ms_j);
    h_Jpz = (REAL*) Malloc(ms_j);
    h_Jmx = (REAL*) Malloc(ms_j);
    h_Jmy = (REAL*) Malloc(ms_j);
    h_Jmz = (REAL*) Malloc(ms_j);
	
    h_Energy = (REAL*) Malloc(ms_s);
	
    h_sdbL = findlog2(h_side);
    h_sdbLs = findlog2(h_side*h_side);
	
    h_sdbLm1 = h_sdbL-1;  //log_2(L) -1 used to divide by L/2
    h_sdbLsm1 = h_sdbLs-1;  //log_2(L*L) -1 used to divide by (L*L)/2
    h_maskL = (h_side>>1)-1;
    h_maskL2 = ((h_side*h_side)>>1)-1;
    h_maskL3 = (V(h_side)>>1)-1;
    h_Ldb2 = h_side>>1;
    h_Lsdb2 = (h_side*h_side)>>1;
	
    h_nspin=h_side*h_side*z_side/2;
	
    srand48(seed);
	
}


//-------------------------------------------------------------------------------------------------------
// CopyConstOnDevice
//-------------------------------------------------------------------------------------------------------
void CopyConstOnDevice(void){
	
    int tmp[2];
	
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_side",&h_side,sizeof(int),0,cudaMemcpyHostToDevice) );
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_nspin",&h_nspin,sizeof(int),0,cudaMemcpyHostToDevice) );
	
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_sdbL",&h_sdbL,sizeof(int),0,cudaMemcpyHostToDevice) );
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_sdbLs",&h_sdbLs,sizeof(int),0,cudaMemcpyHostToDevice) );
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_sdbLm1",&h_sdbLm1,sizeof(int),0,cudaMemcpyHostToDevice) );
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_sdbLsm1",&h_sdbLsm1,sizeof(int),0,cudaMemcpyHostToDevice) );
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_maskL",&h_maskL,sizeof(int),0,cudaMemcpyHostToDevice) );
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_maskL2",&h_maskL2,sizeof(int),0,cudaMemcpyHostToDevice) );
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_maskL3",&h_maskL3,sizeof(int),0,cudaMemcpyHostToDevice) );
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_Ldb2",&h_Ldb2,sizeof(int),0,cudaMemcpyHostToDevice) );
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_Lsdb2",&h_Lsdb2,sizeof(int),0,cudaMemcpyHostToDevice) );
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_h",&h_h,sizeof(REAL),0,cudaMemcpyHostToDevice) );
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_Beta",&h_Beta,sizeof(REAL),0,cudaMemcpyHostToDevice) );
	
    tmp[0]=0;
    tmp[1]=-1;
	
    MY_CUDA_CHECK( cudaMemcpyToSymbol("jindred",tmp,sizeof(tmp),0,cudaMemcpyHostToDevice) );
	
    tmp[0]=-1;
    tmp[1]=0;
	
    MY_CUDA_CHECK( cudaMemcpyToSymbol("jindblue",tmp,sizeof(tmp),0,cudaMemcpyHostToDevice) );
	
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_const_one",&one,sizeof(REAL),0,cudaMemcpyHostToDevice) );
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_const_two",&two,sizeof(REAL),0,cudaMemcpyHostToDevice) );
    MY_CUDA_CHECK( cudaMemcpyToSymbol("d_const_zerotwentyfive",&zerotwentyfive,sizeof(REAL),0,cudaMemcpyHostToDevice) );
	
	
	
}


//-------------------------------------------------------------------------------------------------------
// findlog2
//-------------------------------------------------------------------------------------------------------
int findlog2(int n){
	
    int i=0;
	
    if(n<1 || n&1) { writelog(TRUE,FINDLOG_RC,"Invalid value %d. Must be > 0 and even\n",n); }
	
    while(n>0) {
        n>>=1;
        if(n>0) i++;
    }
	
    return i;
}


//-------------------------------------------------------------------------------------------------------
// nextPow2
//-------------------------------------------------------------------------------------------------------
unsigned int nextPow2( unsigned int x ) {
    --x;
    x |= x >> 1;
    x |= x >> 2;
    x |= x >> 4;
    x |= x >> 8;
    x |= x >> 16;
    return ++x;
}


//-------------------------------------------------------------------------------------------------------
// FreeDevice
//-------------------------------------------------------------------------------------------------------
void FreeDevice(void){
	
    cudaFree(d_Ra);
    cudaFree(d_Rb);
    cudaFree(d_Rc);
    cudaFree(d_Ba);
    cudaFree(d_Bb);
    cudaFree(d_Bc);
	
    cudaFree(d_Jpx);
    cudaFree(d_Jpy);
    cudaFree(d_Jpz);
    cudaFree(d_Jmx);
    cudaFree(d_Jmy);
    cudaFree(d_Jmz);
	
    cudaFree(d_Energy);
	
}


//-------------------------------------------------------------------------------------------------------
// FreeHost
//-------------------------------------------------------------------------------------------------------
void FreeHost(void){
	
    Free(h_Ra);
    Free(h_Rb);
    Free(h_Rc);
    Free(h_Ba);
    Free(h_Bb);
    Free(h_Bc);
	
    Free(h_Jpx);
    Free(h_Jpy);
    Free(h_Jpz);
    Free(h_Jmx);
    Free(h_Jmy);
    Free(h_Jmz);
	
    Free(h_Energy);
	
}


//-------------------------------------------------------------------------------------------------------
// Dump_Spin
//-------------------------------------------------------------------------------------------------------
void Dump_Spin(int count){
	
    FILE* tempfile;
    char namefile[MAXSTRLEN];
	
    return;
    memset(namefile,'\0',MAXSTRLEN); // NON CONTROLLO SE TORNA BENE O MENO
	
    if(count>=0){
		
        snprintf(namefile,sizeof(namefile),"spin_dump_%d_%d",count,mpiid);
		
        MY_CUDA_CHECK( cudaMemcpy( h_Ra, d_Ra+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
        MY_CUDA_CHECK( cudaMemcpy( h_Rb, d_Rb+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
        MY_CUDA_CHECK( cudaMemcpy( h_Rc, d_Rc+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
		
        MY_CUDA_CHECK( cudaMemcpy( h_Ba, d_Ba+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
        MY_CUDA_CHECK( cudaMemcpy( h_Bb, d_Bb+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
        MY_CUDA_CHECK( cudaMemcpy( h_Bc, d_Bc+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
		
    } else {
        snprintf(namefile,sizeof(namefile),"spin_dump_start_%d",mpiid);
    }
	
    tempfile = Fopen(namefile,"w");
	
    for(int i = 0; i<ns_j; i++){
        if(sizeof(REAL)==sizeof(float)) {
            fprintf(tempfile,"%f %f %f\n%f %f %f\n",
                    h_Ra[i],h_Rb[i],h_Rc[i],h_Ba[i],h_Bb[i],h_Bc[i]);
        } else {
            fprintf(tempfile,"%lf %lf %lf\n%lf %lf %lf\n",
                    h_Ra[i],h_Rb[i],h_Rc[i],h_Ba[i],h_Bb[i],h_Bc[i]);
        }
    }
	
    fflush(tempfile);
	
    fclose(tempfile);
	
}
//-------------------------------------------------------------------------------------------------------
// Dump_J
//-------------------------------------------------------------------------------------------------------
void Dump_J(void){
	
    FILE* tempfile;
	
    tempfile = Fopen("j_dump_start","w");
	
    if(sizeof(REAL)==sizeof(float)) {
		
        for(int i = 0; i<ns_j; i++) {
            fprintf(tempfile,"%f %f %f\n",h_Jpx[i],h_Jpy[i],h_Jpz[i]);
            fprintf(tempfile,"%f %f %f\n",h_Jmx[i],h_Jmy[i],h_Jmz[i]);
        }
		
    }else{
		
        for(int i = 0; i<ns_j; i++) {
            fprintf(tempfile,"%lf %lf %lf\n",h_Jpx[i],h_Jpy[i],h_Jpz[i]);
            fprintf(tempfile,"%lf %lf %lf\n",h_Jmx[i],h_Jmy[i],h_Jmz[i]);
        }
		
    }
	
	
    fflush(tempfile);
	
    fclose(tempfile);
	
}


//-------------------------------------------------------------------------------------------------------
// Initialize_Spin
//-------------------------------------------------------------------------------------------------------
void Initialize_Spin(REAL* sa_a, REAL* sa_b, REAL* sa_c){
	
    REAL a, b, c;
    REAL invnorma, norma;
	
    for(int i = 0; i<h_nspin; i++){
		
        norma=two;
		
        while (norma>one){
			
            a=(RANDOM01)*two-one;
            b=(RANDOM01)*two-one;
            c=(RANDOM01)*two-one;
			
            norma=a*a+b*b+c*c;
			
        }
		
        invnorma=one/SQRT(norma);
		
        sa_a[i]=a*invnorma;
        sa_b[i]=b*invnorma;
        sa_c[i]=c*invnorma;
		
    }
	
}


//-------------------------------------------------------------------------------------------------------
// Read_Spin
//-------------------------------------------------------------------------------------------------------
void Read_Spin(char *inputfile,
               REAL* sr_a, REAL* sr_b, REAL* sr_c,
               REAL* sb_a, REAL* sb_b, REAL* sb_c) {
	
    REAL a, b, c;
    int j=0, k=0;
    FILE *fpi;
    fpi=Fopen(inputfile,"r");
    fprintf(stderr,"Reading spin from file %s\n",inputfile);
	
    for(int i = 0; i< V(h_side); i++){
        if(sizeof(REAL)==sizeof(float)) {
            fscanf(fpi,"%f %f %f",&a,&b,&c);
        } else {
            fscanf(fpi,"%lf %lf %lf",&a,&b,&c);
        }
        if(i&1) {
            sb_a[j]=a;
            sb_b[j]=b;
            sb_c[j]=c;
            j++;
        } else {
            sr_a[k]=a;
            sr_b[k]=b;
            sr_c[k]=c;
            k++;
        }
    }
    fclose(fpi);
	
}


//-------------------------------------------------------------------------------------------------------
// CopyDataOnDevice
//-------------------------------------------------------------------------------------------------------
void CopyDataOnDevice(void){
	
    MY_CUDA_CHECK( cudaMalloc( (void **) &d_Ra, ms_s_mpi ) );
    MY_CUDA_CHECK( cudaMalloc( (void **) &d_Rb, ms_s_mpi ) );
    MY_CUDA_CHECK( cudaMalloc( (void **) &d_Rc, ms_s_mpi ) );
	
    MY_CUDA_CHECK( cudaMalloc( (void **) &d_Ba, ms_s_mpi ) );
    MY_CUDA_CHECK( cudaMalloc( (void **) &d_Bb, ms_s_mpi ) );
    MY_CUDA_CHECK( cudaMalloc( (void **) &d_Bc, ms_s_mpi ) );
	
    MY_CUDA_CHECK( cudaMemcpy( d_Ra+h_Lsdb2, h_Ra, ms_s, cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Rb+h_Lsdb2, h_Rb, ms_s, cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Rc+h_Lsdb2, h_Rc, ms_s, cudaMemcpyHostToDevice ) );
	
    MY_CUDA_CHECK( cudaMemcpy( d_Ba+h_Lsdb2, h_Ba, ms_s, cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Bb+h_Lsdb2, h_Bb, ms_s, cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Bc+h_Lsdb2, h_Bc, ms_s, cudaMemcpyHostToDevice ) );
	
    MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jpx, ms_j_mpi ) );
    MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jpy, ms_j_mpi ) );
    MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jpz, ms_j_mpi ) );
	
    MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jmx, ms_j_mpi ) );
    MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jmy, ms_j_mpi ) );
    MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jmz, ms_j_mpi ) );
	
    MY_CUDA_CHECK( cudaMemcpy( d_Jpx+h_Lsdb2, h_Jpx, ms_j, cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Jpy+h_Lsdb2, h_Jpy, ms_j, cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Jpz+h_Lsdb2, h_Jpz, ms_j, cudaMemcpyHostToDevice ) );
	
    MY_CUDA_CHECK( cudaMemcpy( d_Jmx+h_Lsdb2, h_Jmx, ms_j, cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Jmy+h_Lsdb2, h_Jmy, ms_j, cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Jmz+h_Lsdb2, h_Jmz, ms_j, cudaMemcpyHostToDevice ) );
	
    MY_CUDA_CHECK( cudaMalloc( (void **) &d_Energy, ms_s ) );
    MY_CUDA_CHECK( cudaMemset( d_Energy, 0, ms_s ) );
}


//-------------------------------------------------------------------------------------------------------
// Initialize_J
//-------------------------------------------------------------------------------------------------------
void Initialize_J(void){
	
    int j;
	
    switch (init_type_j){
		
    case 0:
        for (j=0;j<ns_j;j++) {
            h_Jpx[j] = ((RANDOM01)>zerofive) ? -1 : 1;
            h_Jpy[j] = ((RANDOM01)>zerofive) ? -1 : 1;
            h_Jpz[j] = ((RANDOM01)>zerofive) ? -1 : 1;
            h_Jmx[j] = ((RANDOM01)>zerofive) ? -1 : 1;
            h_Jmy[j] = ((RANDOM01)>zerofive) ? -1 : 1;
            h_Jmz[j] = ((RANDOM01)>zerofive) ? -1 : 1;
        }
		
        break;
    case 1:
        for (j=0;j<ns_j;j++) {
            h_Jpx[j] = h_Jpy[j] = h_Jpz[j] =
                h_Jmx[j] = h_Jmy[j] = h_Jmz[j] = one;
        }
        break;
    case -1:
        for (j=0;j<ns_j;j++) {
            h_Jpx[j] = h_Jpy[j] = h_Jpz[j] =
                h_Jmx[j] = h_Jmy[j] = h_Jmz[j] = -one;
        }
		
        break;
    case 2: //Caso Gaussiano
        for (j=0;j<ns_j;j++) {
            h_Jpx[j] = rand_gauss();
            h_Jpy[j] = rand_gauss();
            h_Jpz[j] = rand_gauss();
            h_Jmx[j] = rand_gauss();
            h_Jmy[j] = rand_gauss();
            h_Jmz[j] = rand_gauss();
        }
        break;
    default: {
        writelog(TRUE,APPLICATION_RC,"Invalid initialization for J: %d", init_type_j);
        break;
    }
    }
	
}


//-------------------------------------------------------------------------------------------------------
// rand_gauss
//-------------------------------------------------------------------------------------------------------
REAL rand_gauss(void){
	
    REAL r,cosr;
    REAL a=one-RANDOM01;
    if(a<EPSILON) { a=EPSILON; }
	
    r=SQRT(-2.*LOG(a));
	
    cosr=COS((2*M_PI)*(RANDOM01));
	
    return r*cosr;
	
}


//-------------------------------------------------------------------------------------------------------
// Parameters_Reading
//-------------------------------------------------------------------------------------------------------
void Parameters_Reading(char *inputfile){
	
    dictionary *ini;
    char key[MAXSTRLEN];
	
    if(inputfile==NULL) { writelog(TRUE,APPLICATION_RC,"No inputfile"); }
	
    ini = iniparser_load(inputfile);
	
    if(ini==NULL) { writelog(TRUE,APPLICATION_RC,"Cannot parse file: %s\n", inputfile); }
	
    READINTFI(h_side,"Side");
    READINTFI(init_type_j,"Initialization type for J");
    READINTFI(nsteps,"Number of steps");
    READINTFI(overrelaxiter,"Number of overrelaxation iterations");
    READINTFI(outfreq,"Output frequency");
    READINTFI(dumpfreq,"Dump configuration frequency");
    READREALFI(h_h,"Magnetic field");
    READREALFI(h_Beta,"Beta (1/T)");
    READINTFI(seed,"Random seed");
	
    DEVICEREADINTFI(numBlocks,"Number of blocks");
    DEVICEREADINTFI(numThreadsPerBlock,"Number of threads for block");
    DEVICEREADINTFI(whichgpu,"Which GPU");
	
    DEVICEREADINTFI(RandNumBlocks,"RandGene num of blocks");
    DEVICEREADINTFI(RandNumThreadsPerBlock,"RandGene num of threads for block");
	
}


//-------------------------------------------------------------------------------------------------------
// Calculate_Sizes
//-------------------------------------------------------------------------------------------------------
void Calculate_Sizes(void){
	
    ns_j = h_side*h_side*z_side/2;
	
    ms_s = (h_side*h_side*z_side/2)*sizeof(REAL);
    ms_s_mpi = ((h_side*h_side*z_side/2)+(h_side*h_side))*sizeof(REAL);
	
    ms_j = ns_j*sizeof(REAL);
    ms_j_mpi = ((h_side*h_side*z_side/2)+(h_side*h_side))*sizeof(REAL);
	
}


//-------------------------------------------------------------------------------------------------------
// getNumBlocksAndThreads
//-------------------------------------------------------------------------------------------------------
void getNumBlocksAndThreads(int n, int maxBlocks, int maxThreads, int &blocks, int &threads) {
	
    if (n == 1){
        threads = 1;
        blocks = 1;
    } else{
        threads = (n < maxThreads*2) ? nextPow2(n / 2) : maxThreads;
        blocks = max(1, n / (threads * 2));
    }
	
    blocks = min(maxBlocks, blocks);
}


//-------------------------------------------------------------------------------------------------------
// main
//-------------------------------------------------------------------------------------------------------
int main(int argc, char **argv){
	
    char *inputfile = NULL;
    char *po;
    int istep;
    int nblocks_red, nthreads_red;
    int cudadevice, ngpu;
    double ttimeover=0., ttimehb=0.;
    REAL normene;
    int j;
    char *readspin=NULL;
#if !defined(HOSTREDU)
    unsigned int retCnt;
#endif
    int myid=0, nprocs=1;
    int lpoff;
    struct mpiexch *dataexch;
    int numBlocksBnd;
	
    cudaStream_t *stream;
    int nstreams=2;
	
    TIMER_DEF;
	
    for(int i = 1; i < argc; i++) {
        po = argv[i];
        if (*po++ == '-') {
            switch (*po++) {
            case '?':
                SKIPBLANK
                    Usage(argv[0]);
                exit(OK);
                break;
            case 'v':
                SKIPBLANK
                    verbose=TRUE;
                break;
            case 'i':
                SKIPBLANK
                    inputfile=Strdup(po);
                break;
            case 'o':
                SKIPBLANK
                    onlyOverrelaxation=true;
                break;
            case 'r':
                SKIPBLANK
                    readspin=Strdup(po);
                break;
				
            case 'b':
                SKIPBLANK
                    onlyHeatBath=true;
                break;
            default:
                Usage(argv[0]);
                exit(OK);
                break;
            }
        }
    }
    if(inputfile==NULL) {
        Usage(argv[0]);
        exit(OK);
    }
	
    StartMpi(&myid, &nprocs,&argc,&argv);
    mpiid=myid;
	
#if defined(DEBUG)
    cudaPrintfInit();
#endif

    cudaGetDeviceCount(&ngpu);
	
    if(ngpu<1) {
        writelog(TRUE,APPLICATION_RC,"No gpu available! %d\n",ngpu);
    }
	
    Parameters_Reading(inputfile);
	
    seed=(seed+myid*nprocs)%MAXSEED;
	
	
    z_side=h_side/nprocs;
	
    numBlocksBnd=h_Lsdb2/numThreadsPerBlock+(h_Lsdb2%numThreadsPerBlock?1:0);
	
    if((z_side*nprocs)<h_side) {
        writelog(TRUE,APPLICATION_RC,"Invalid number of processors (%d) or system size (%d). Number of processors must be an exact divider of the system size", nprocs, h_side);
    }
	
    if(whichgpu>(ngpu-1)) {
        writelog(TRUE,APPLICATION_RC, 
                 "invalid number of gpu: %d must be > 0 and <= %d\n",whichgpu, ngpu);
    }
	
	
    if(nprocs>1) {
        if(myid%2) { /* Odd task go on GPU 1 */
            whichgpu=1;
        } else {
            whichgpu=0;
        }
    }

    if(whichgpu>=0) {
        printf("There are %d GPUs!, forcing usage of n. <%d>\n",ngpu,whichgpu);
        MY_CUDA_CHECK(cudaSetDevice(whichgpu));
    }

    dataexch=SetUpMpi(myid,nprocs,h_side);
	
    if (h_side<=0 || !pow_of_two(h_side)) {
        writelog(TRUE,APPLICATION_RC,"Invalid side: must be > 0 and power of two\n");
    }
	
    Calculate_Sizes();
    Init();
	
    lpoff=(z_side-1)*h_Lsdb2;
	
    if(readspin==NULL) {
		
        Initialize_Spin(h_Ra,h_Rb,h_Rc);
        Initialize_Spin(h_Ba,h_Bb,h_Bc);
        Dump_Spin(-1);
		
    } else {
        Read_Spin(readspin,h_Ra,h_Rb,h_Rc, h_Ba,h_Bb,h_Bc);
    }
	
    Initialize_J();
	
	
    // Dump_J();
	
    initRand(V(h_side));
	
    CopyConstOnDevice();
	
    CopyDataOnDevice();
	
    getNumBlocksAndThreads(h_nspin, numBlocks, numThreadsPerBlock,nblocks_red, nthreads_red);
	
    normene=one/(three*V(h_side));
	
    h_totene=zeroEne;
	
    if(myid==0 && nprocs==1) {
        ComputeHEnergy(h_Ra, h_Rb, h_Rc, h_Ba, h_Bb, h_Bc,h_Jpx, h_Jpy, h_Jpz, h_Jmx, h_Jmy, h_Jmz, &h_totene,TRUE);
        h_totene=SumEnergy(h_totene, myid);
        printf("Initial Energy per spin on host: %f\n", normene*h_totene);
    }
	
    memcpy(dataexch->downsend,h_Ra,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->downsend+h_Lsdb2,h_Rb,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->downsend+(2*h_Lsdb2),h_Rc,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->downsend+(3*h_Lsdb2),h_Ba,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->downsend+(4*h_Lsdb2),h_Bb,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->downsend+(5*h_Lsdb2),h_Bc,h_Lsdb2*sizeof(REAL));
	
    memcpy(dataexch->upsend,h_Ra+lpoff,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->upsend+h_Lsdb2,h_Rb+lpoff,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->upsend+(2*h_Lsdb2),h_Rc+lpoff,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->upsend+(3*h_Lsdb2),h_Ba+lpoff,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->upsend+(4*h_Lsdb2),h_Bb+lpoff,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->upsend+(5*h_Lsdb2),h_Bc+lpoff,h_Lsdb2*sizeof(REAL));
	
	
    dataexch->nbyte=2*h_Lsdb2*sizeof(REAL)*NCOMP;
    dataexch->post=FALSE;
	
    ExchMpi(dataexch);
	
    MY_CUDA_CHECK( cudaMemcpy( d_Ra, dataexch->downrecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Rb, dataexch->downrecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Rc, dataexch->downrecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Ba, dataexch->downrecv+(3*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Bb, dataexch->downrecv+(4*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Bc, dataexch->downrecv+(5*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	
    MY_CUDA_CHECK( cudaMemcpy( d_Ra+h_nspin+h_Lsdb2, dataexch->uprecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Rb+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Rc+h_nspin+h_Lsdb2, dataexch->uprecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Ba+h_nspin+h_Lsdb2, dataexch->uprecv+(3*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Bb+h_nspin+h_Lsdb2, dataexch->uprecv+(4*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Bc+h_nspin+h_Lsdb2, dataexch->uprecv+(5*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	
    memcpy(dataexch->downsend,h_Jpx,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->downsend+h_Lsdb2,h_Jpy,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->downsend+(2*h_Lsdb2),h_Jpz,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->downsend+(3*h_Lsdb2),h_Jmx,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->downsend+(4*h_Lsdb2),h_Jmy,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->downsend+(5*h_Lsdb2),h_Jmz,h_Lsdb2*sizeof(REAL));
	
    memcpy(dataexch->upsend,h_Jpx+lpoff,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->upsend+h_Lsdb2,h_Jpy+lpoff,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->upsend+(2*h_Lsdb2),h_Jpz+lpoff,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->upsend+(3*h_Lsdb2),h_Jmx+lpoff,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->upsend+(4*h_Lsdb2),h_Jmy+lpoff,h_Lsdb2*sizeof(REAL));
    memcpy(dataexch->upsend+(5*h_Lsdb2),h_Jmz+lpoff,h_Lsdb2*sizeof(REAL));
	
    dataexch->post=FALSE;
	
    ExchMpi(dataexch);
	
    MY_CUDA_CHECK( cudaMemcpy( d_Jpx, dataexch->downrecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Jpy, dataexch->downrecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Jpz, dataexch->downrecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Jmx, dataexch->downrecv+(3*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Jmy, dataexch->downrecv+(4*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Jmz, dataexch->downrecv+(5*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	
    MY_CUDA_CHECK( cudaMemcpy( d_Jpx+h_nspin+h_Lsdb2, dataexch->uprecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Jpy+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Jpz+h_nspin+h_Lsdb2, dataexch->uprecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Jmx+h_nspin+h_Lsdb2, dataexch->uprecv+(3*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Jmy+h_nspin+h_Lsdb2, dataexch->uprecv+(4*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
    MY_CUDA_CHECK( cudaMemcpy( d_Jmz+h_nspin+h_Lsdb2, dataexch->uprecv+(5*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
	
    dataexch->post=FALSE;
	
    dataexch->nbyte=h_Lsdb2*sizeof(REAL)*NCOMP;
	
    MY_CUDA_CHECK( cudaMalloc((void **)&r_odata,numBlocks*sizeof(REAL)) );
	
    dim3 dimGrid(numBlocks,1,1);
    dim3 dimGridBnd(numBlocksBnd,1,1);
    dim3 dimBlock(numThreadsPerBlock,1,1);
	
    ComputeEnergy<<< dimGrid,dimBlock >>>(d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                          d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                          d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
                                          d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
                                          d_Energy, TRUE);
										  
#if defined(HOSTREDU)
    MY_CUDA_CHECK( cudaMemcpy(h_Energy, d_Energy, h_nspin*sizeof(REAL),cudaMemcpyDeviceToHost) );
	
    h_totene=zeroEne;
	
    for(int i=0; i<h_nspin; i++) { h_totene+=h_Energy[i]; }
	
#else

    MY_CUDA_CHECK( cudaThreadSynchronize() );
    retCnt=0;
    cudaMemcpyToSymbol("retirementCount", &retCnt, sizeof(unsigned int), 0, cudaMemcpyHostToDevice);
    reduceSinglePass(h_nspin, nthreads_red, nblocks_red, d_Energy, r_odata);
    MY_CUDA_CHECK( cudaThreadSynchronize() );
    MY_CUDA_CHECK( cudaMemcpy(&h_totene, r_odata, sizeof(REAL),cudaMemcpyDeviceToHost) );
	
#endif

    h_totene=SumEnergy(h_totene, myid);
    if(myid==0) {
        printf("Initial Energy per spin: %f\n", normene*h_totene);
    }
	
#if defined(FERMI)
    cudaFuncSetCacheConfig(OverrelaxRedBnd, cudaFuncCachePreferL1 );
    cudaFuncSetCacheConfig(OverrelaxRedBulk, cudaFuncCachePreferL1 );
    cudaFuncSetCacheConfig(OverrelaxBlueBnd, cudaFuncCachePreferL1 );
    cudaFuncSetCacheConfig(OverrelaxBlueBulk, cudaFuncCachePreferL1 );
#endif
    // allocate and initialize an array of stream handles
    stream = (cudaStream_t*) Malloc(nstreams * sizeof(cudaStream_t));
	
    for(int i = 0; i < nstreams; i++) {
        MY_CUDA_CHECK( cudaStreamCreate(&(stream[i])) );
    }
		
    for(istep=0; istep<nsteps; istep++) {
		
        if(!onlyHeatBath){
			
            TIMER_START;
			
            for(j=0; j<overrelaxiter; j++) {
				
#if defined(SINGLEOVERKERNEL)

                Overrelaxation<<< dimGrid,dimBlock>>>(d_Ra, d_Rb, d_Rc,
                                                      d_Ba, d_Bb, d_Bc,
                                                      d_Jpx, d_Jpy, d_Jpz,
                                                      d_Jmx, d_Jmy, d_Jmz, TRUE);
				
                MY_CUDA_CHECK( cudaThreadSynchronize() );
				
                Overrelaxation<<< dimGrid,dimBlock>>>(d_Ba, d_Bb, d_Bc,
                                                      d_Ra, d_Rb, d_Rc,
                                                      d_Jmx, d_Jmy, d_Jmz,
                                                      d_Jpx, d_Jpy, d_Jpz, FALSE);
				
                MY_CUDA_CHECK( cudaThreadSynchronize() );
				
#else

                OverrelaxRedBnd<<< dimGridBnd,dimBlock, 0, stream[1] >>>(d_Ra+h_Lsdb2,
                                                                         d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                                                         d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                                                         d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
                                                                         d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2);
                OverrelaxRedBulk<<< dimGrid,dimBlock, 0, stream[0] >>>(d_Ra+h_Lsdb2,
                                                                       d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                                                       d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                                                       d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
                                                                       d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2);
																	   
                MY_CUDA_CHECK(cudaStreamSynchronize(stream[0]));
				
                MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend, d_Ra+h_Lsdb2,h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+h_Lsdb2, d_Rb+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+(2*h_Lsdb2), d_Rc+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
                MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend, d_Ra+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+h_Lsdb2, d_Rb+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+(2*h_Lsdb2), d_Rc+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
                MY_CUDA_CHECK(cudaStreamSynchronize(stream[1]));
				
                ExchMpi(dataexch);
				
                MY_CUDA_CHECK( cudaMemcpyAsync( d_Ra, dataexch->downrecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( d_Rb, dataexch->downrecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( d_Rc, dataexch->downrecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				
                MY_CUDA_CHECK( cudaMemcpyAsync( d_Ra+h_nspin+h_Lsdb2, dataexch->uprecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( d_Rb+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( d_Rc+h_nspin+h_Lsdb2,dataexch->uprecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				
                MY_CUDA_CHECK( cudaThreadSynchronize() );
				
                OverrelaxBlueBnd<<< dimGridBnd,dimBlock, 0, stream[1]>>>(d_Ba+h_Lsdb2,
                                                                         d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                                                         d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                                                         d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
                                                                         d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2);
                OverrelaxBlueBulk<<< dimGrid,dimBlock, 0, stream[0]>>>(d_Ba+h_Lsdb2,
                                                                       d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                                                       d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                                                       d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
                                                                       d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2);
																	   
                MY_CUDA_CHECK(cudaStreamSynchronize(stream[0]));
				
                MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend, d_Ba+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+h_Lsdb2, d_Bb+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+(2*h_Lsdb2), d_Bc+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
                MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend, d_Ba+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+h_Lsdb2, d_Bb+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+(2*h_Lsdb2), d_Bc+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
                MY_CUDA_CHECK(cudaStreamSynchronize(stream[1]));
				
                ExchMpi(dataexch);
				
                MY_CUDA_CHECK( cudaMemcpyAsync( d_Ba, dataexch->downrecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( d_Bb, dataexch->downrecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( d_Bc, dataexch->downrecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				
                MY_CUDA_CHECK( cudaMemcpyAsync( d_Ba+h_nspin+h_Lsdb2, dataexch->uprecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( d_Bb+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
                MY_CUDA_CHECK( cudaMemcpyAsync( d_Bc+h_nspin+h_Lsdb2,dataexch->uprecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
                MY_CUDA_CHECK( cudaThreadSynchronize() );
				
#endif
            }
			
            TIMER_STOP;
            ttimeover+=TIMER_ELAPSED;
        }
		
        if(!onlyOverrelaxation) {
			
            TIMER_START;
			
            /* Here it goes the heat-bath */
            /* first of all, random number generation */
            /* then the heat-bath itself */
            ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ONEONE,d_cRand);
            ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ZEROFIVE,d_aRand);
            ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ZEROFIVE,d_bRand);
            ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ZEROONE,d_Rand);
			      
            if(istep==0&&FALSE){
                printRand(randNum,d_aRand,"daRand3");
                printRand(randNum,d_bRand,"dbRand3");
                printRand(randNum,d_cRand,"dcRand3");
                printRand(randNum,d_Rand,"dRand3");
                exit(1);
            }
			
			
#if defined(SINGLEHBKERNEL)

            HeatBath<<< dimGrid,dimBlock >>>(d_Ra, d_Rb, d_Rc,
                                             d_Ba, d_Bb, d_Bc,
                                             d_Jpx, d_Jpy, d_Jpz,
                                             d_Jmx, d_Jmy, d_Jmz,
                                             d_aRand,d_bRand,d_cRand,
                                             d_Rand,TRUE);
			
            MY_CUDA_CHECK( cudaThreadSynchronize() );
			
            HeatBath<<< dimGrid,dimBlock >>>(d_Ba, d_Bb, d_Bc,
                                             d_Ra, d_Rb, d_Rc,
                                             d_Jmx, d_Jmy, d_Jmz,
                                             d_Jpx, d_Jpy, d_Jpz,
                                             d_aRand+h_nspin,d_bRand+h_nspin,
                                             d_cRand+h_nspin,
                                             d_Rand+h_nspin,FALSE);
			
            MY_CUDA_CHECK( cudaThreadSynchronize() );
      
#else


            HeatBathRedBnd<<< dimGridBnd,dimBlock, 0, stream[1] >>>(d_Ra+h_Lsdb2,
                                                                    d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                                                    d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                                                    d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
                                                                    d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
                                                                    d_aRand,d_bRand,d_cRand,d_Rand);
							
            HeatBathRedBulk<<< dimGrid,dimBlock, 0, stream[0] >>>(d_Ra+h_Lsdb2,
                                                                  d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                                                  d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                                                  d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
                                                                  d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
                                                                  d_aRand,d_bRand,d_cRand,d_Rand);
																	   
            MY_CUDA_CHECK(cudaStreamSynchronize(stream[0]));

            MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend, d_Ra+h_Lsdb2,h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+h_Lsdb2, d_Rb+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+(2*h_Lsdb2), d_Rc+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
            MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend, d_Ra+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+h_Lsdb2, d_Rb+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+(2*h_Lsdb2), d_Rc+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
            MY_CUDA_CHECK(cudaStreamSynchronize(stream[1]));
				
            ExchMpi(dataexch);
				
            MY_CUDA_CHECK( cudaMemcpyAsync( d_Ra, dataexch->downrecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( d_Rb, dataexch->downrecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( d_Rc, dataexch->downrecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				
            MY_CUDA_CHECK( cudaMemcpyAsync( d_Ra+h_nspin+h_Lsdb2, dataexch->uprecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( d_Rb+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( d_Rc+h_nspin+h_Lsdb2,dataexch->uprecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				
            MY_CUDA_CHECK( cudaThreadSynchronize() );
				
            HeatBathBlueBnd<<< dimGridBnd,dimBlock, 0, stream[1]>>>(d_Ba+h_Lsdb2,
                                                                    d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                                                    d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                                                    d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
                                                                    d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
                                                                    d_aRand,d_bRand,d_cRand,d_Rand);

            HeatBathBlueBulk<<< dimGrid,dimBlock, 0, stream[0]>>>(d_Ba+h_Lsdb2,
                                                                  d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                                                  d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                                                  d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
                                                                  d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
                                                                  d_aRand,d_bRand,d_cRand,d_Rand);
																	   
            MY_CUDA_CHECK(cudaStreamSynchronize(stream[0]));

            MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend, d_Ba+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+h_Lsdb2, d_Bb+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->downsend+(2*h_Lsdb2), d_Bc+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
            MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend, d_Ba+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+h_Lsdb2, d_Bb+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( dataexch->upsend+(2*h_Lsdb2), d_Bc+h_nspin, h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost, stream[1] ) );
				
            MY_CUDA_CHECK(cudaStreamSynchronize(stream[1]));
				
            ExchMpi(dataexch);
				
            MY_CUDA_CHECK( cudaMemcpyAsync( d_Ba, dataexch->downrecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( d_Bb, dataexch->downrecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( d_Bc, dataexch->downrecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
				
            MY_CUDA_CHECK( cudaMemcpyAsync( d_Ba+h_nspin+h_Lsdb2, dataexch->uprecv, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( d_Bb+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2, h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );
            MY_CUDA_CHECK( cudaMemcpyAsync( d_Bc+h_nspin+h_Lsdb2,dataexch->uprecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice, stream[1] ) );

            MY_CUDA_CHECK( cudaThreadSynchronize() );




            /*		HeatBathRed<<< dimGrid,dimBlock >>>(d_Ra, d_Rb, d_Rc,
                    d_Ba, d_Bb, d_Bc,
                    d_Jpx, d_Jpy, d_Jpz,
                    d_Jmx, d_Jmy, d_Jmz,
                    d_aRand,d_bRand,d_cRand,
                    d_Rand);
			
                    MY_CUDA_CHECK( cudaThreadSynchronize() );
			
                    HeatBathBlue<<< dimGrid,dimBlock >>>(d_Ba, d_Bb, d_Bc,
                    d_Ra, d_Rb, d_Rc,
                    d_Jmx, d_Jmy, d_Jmz,
                    d_Jpx, d_Jpy, d_Jpz,
                    d_aRand+h_nspin,d_bRand+h_nspin,
                    d_cRand+h_nspin,
                    d_Rand+h_nspin);
			
                    MY_CUDA_CHECK( cudaThreadSynchronize() );*/
			
#endif
			
            TIMER_STOP;
            ttimehb+=TIMER_ELAPSED;
			
        }
		
        if(istep && (istep%outfreq)==0) {
            ComputeEnergy<<< dimGrid,dimBlock >>>(d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                                  d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                                  d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
                                                  d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
                                                  d_Energy, TRUE);
			
            MY_CUDA_CHECK( cudaThreadSynchronize() );
			
			
#if defined(HOSTREDU)
			
            MY_CUDA_CHECK( cudaMemcpy(h_Energy, d_Energy,
                                      h_nspin*sizeof(REAL),
                                      cudaMemcpyDeviceToHost) );
            h_totene=zeroEne;
			
            for(int i=0; i<h_nspin; i++) h_totene+=h_Energy[i];
			
#else
            retCnt=0;
            cudaMemcpyToSymbol("retirementCount", &retCnt, sizeof(unsigned int), 0, cudaMemcpyHostToDevice);
            reduceSinglePass(h_nspin, nthreads_red, nblocks_red, d_Energy, r_odata);
            MY_CUDA_CHECK( cudaMemcpy(&h_totene, r_odata, sizeof(REAL),cudaMemcpyDeviceToHost) );
			
#endif
            h_totene=SumEnergy(h_totene, myid);
            if(myid==0) {
                printf("Energy per spin at step %d: %f\n",istep,normene*h_totene);
            }
            //place holder for energy evaluation
        }
		
        if(istep && (istep%dumpfreq)==0) {
            //Dump_Spin(istep);
        }
		
    }
	
    if(myid==0) {
        if(!onlyHeatBath){
            printf("Total time for %d iterations of overrelaxation on a %d^3 cube: %f seconds\n", nsteps*overrelaxiter, h_side, ttimeover*0.000001);
            printf("Time for spin: %f microseconds\n",ttimeover/(overrelaxiter*(double)nsteps*V(h_side)));
        }
		
        if(!onlyOverrelaxation){
            printf("Total time for %d iterations of heat bath on a %d^3 cube: %f seconds\n", nsteps, h_side, ttimehb*0.000001);
            printf("Time for spin: %f microseconds\n",ttimehb/((double)nsteps*V(h_side)));
        }
    }
	
    ComputeEnergy<<< dimGrid,dimBlock >>>(d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                          d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                          d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
                                          d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
                                          d_Energy, TRUE);
	
#if defined(HOSTREDU)
	
    MY_CUDA_CHECK( cudaMemcpy(h_Energy, d_Energy, h_nspin*sizeof(REAL),cudaMemcpyDeviceToHost) );
	
    h_totene=zero;
	
    for(int i=0; i<h_nspin; i++) h_totene+=h_Energy[i];
	
#else
	
    MY_CUDA_CHECK( cudaThreadSynchronize() );
    reduceSinglePass(ns_j, nthreads_red, nblocks_red,d_Energy, r_odata);
    MY_CUDA_CHECK( cudaMemcpy(&h_totene, r_odata, sizeof(REAL),cudaMemcpyDeviceToHost) );
	
#endif

    h_totene=SumEnergy(h_totene, myid);
	
    if(myid==0) {
        printf("Final Energy per spin: %f\n",normene*h_totene);
    }
	
    Dump_Spin(nsteps);
	
    FreeDevice();
	
    FreeHost();
	
    FreeRand();
	
#if defined(DEBUG)
    cudaPrintfDisplay();
    cudaPrintfEnd();
#endif

    StopMpi();
	
    return 0;
	
}


//-------------------------------------------------------------------------------------------------------
// Overrelaxation
//-------------------------------------------------------------------------------------------------------
#if defined(SINGLEOVERKERNEL)
__global__ void Overrelaxation(REAL *a, REAL *b, REAL *c,
                               REAL *an, REAL *bn, REAL *cn,
                               REAL *Jpx, REAL *Jpy, REAL *Jpz,
                               REAL *Jmx, REAL *Jmy, REAL *Jmz,
                               int even) {
	
    unsigned int iy, iz;
    int diff;
    REAL as, bs, cs;
	
    int spx, smx, spy, smy, spz, smz;
	
    REAL factor;
	
    int i;
	
    const unsigned int tid = threadIdx.x;
	
    const unsigned int num_threads = gridDim.x*blockDim.x;
    int istart;
	
    istart=blockDim.x*blockIdx.x+d_Lsdb2; /*start on second plane (by using d_Lsdb2)
                                            since first plane is for the boundary */
	
	
    for(i=istart+tid; i<d_nspin; i+=num_threads) {
		
        iz=i>>d_sdbLsm1;
        iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
        //do the same thing if both $iz and $iy are even or odd
        diff=(((iz&1)^(iy&1))^even);
		
        spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
        smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
        spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
        smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
        spz=(i+d_Lsdb2) & d_maskL3;
        smz=(i-d_Lsdb2) & d_maskL3;
		
        as=an[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
            an[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
            an[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
            an[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
            an[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
            an[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];
		
        bs=bn[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
            bn[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
            bn[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
            bn[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
            bn[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
            bn[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];
		
        cs=cn[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
            cn[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
            cn[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
            cn[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
            cn[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
            cn[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];
		
        cs+=d_h;
		
        //Overrelaxation
		
        factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/(as*as+bs*bs+cs*cs);
		
        a[i]=as*factor-a[i];
        b[i]=bs*factor-b[i];
        c[i]=cs*factor-c[i];
		
    }
	
}

#else

//-------------------------------------------------------------------------------------------------------
// OverrelaxRedBnd
//-------------------------------------------------------------------------------------------------------
__global__ void OverrelaxRedBnd(REAL *a, REAL *b, REAL *c,
                                REAL *an, REAL *bn, REAL *cn,
                                REAL *Jpx, REAL *Jpy, REAL *Jpz,
                                REAL *Jmx, REAL *Jmy, REAL *Jmz) {
	
    unsigned int iy, iz;
    int diff;
    REAL as, bs, cs;
	
    int spx, smx, spy, smy, spz, smz;
	
    REAL factor;
	
    int i;
	
    const unsigned int tid = threadIdx.x;
	
    const unsigned int num_threads = gridDim.x*blockDim.x;
    int istart;
	
    istart=blockDim.x*blockIdx.x; /* does first plane (up to d_Lsdb2) */
	
    for(i=istart+tid; i<d_Lsdb2; i+=num_threads) {
		
        iz=i>>d_sdbLsm1;
        iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
        //do the same thing if both $iz and $iy are even or odd
        diff=(iz&1)^(iy&1);
		
        spx=(i >> (-diff)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (diff-1));
        smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (-diff)) + (i>>(diff-1));
        spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
        smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
        spz=(i+d_Lsdb2) /* & d_maskL3 */;
        smz=(i-d_Lsdb2) /* & d_maskL3 */;
		
        as=an[spx]*Jpx[i]+
            an[spy]*Jpy[i]+
            an[spz]*Jpz[i]+
            an[smx]*Jmx[i]+
            an[smy]*Jmy[i]+
            an[smz]*Jmz[i];
		
        bs=bn[spx]*Jpx[i]+
            bn[spy]*Jpy[i]+
            bn[spz]*Jpz[i]+
            bn[smx]*Jmx[i]+
            bn[smy]*Jmy[i]+
            bn[smz]*Jmz[i];
		
        cs=cn[spx]*Jpx[i]+
            cn[spy]*Jpy[i]+
            cn[spz]*Jpz[i]+
            cn[smx]*Jmx[i]+
            cn[smy]*Jmy[i]+
            cn[smz]*Jmz[i];
		
        cs+=d_h;
		
        //Overrelaxation
		
        factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/(as*as+bs*bs+cs*cs);
		
#if defined(DEBUG)
        printf("Reds [%d]: %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f ",i,a[i],b[i],c[i], as,bs,cs, factor);
#endif
		
        a[i]=as*factor-a[i];
        b[i]=bs*factor-b[i];
        c[i]=cs*factor-c[i];
		
#if defined(DEBUG)
        printf("out: %7.5f %7.5f\n",a[i],b[i],c[i]);
#endif
		
		
    }
    istart=(d_nspin-d_Lsdb2)+blockDim.x*blockIdx.x;
	
    for(i=istart+tid; i<d_nspin; i+=num_threads) {
        iz=i>>d_sdbLsm1;
        iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
        //do the same thing if both $iz and $iy are even or odd
        diff=(iz&1)^(iy&1);
		
        spx=(i >> (-diff)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (diff-1));
        smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (-diff)) + (i>>(diff-1));
        spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
        smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
        spz=(i+d_Lsdb2) /* & d_maskL3 */;
        smz=(i-d_Lsdb2) /* & d_maskL3 */;
		
        as=an[spx]*Jpx[i]+
            an[spy]*Jpy[i]+
            an[spz]*Jpz[i]+
            an[smx]*Jmx[i]+
            an[smy]*Jmy[i]+
            an[smz]*Jmz[i];
		
        bs=bn[spx]*Jpx[i]+
            bn[spy]*Jpy[i]+
            bn[spz]*Jpz[i]+
            bn[smx]*Jmx[i]+
            bn[smy]*Jmy[i]+
            bn[smz]*Jmz[i];
		
        cs=cn[spx]*Jpx[i]+
            cn[spy]*Jpy[i]+
            cn[spz]*Jpz[i]+
            cn[smx]*Jmx[i]+
            cn[smy]*Jmy[i]+
            cn[smz]*Jmz[i];
		
        cs+=d_h;
		
        //Overrelaxation
		
        factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/(as*as+bs*bs+cs*cs);
		
#if defined(DEBUG)
        printf("Reds [%d]: %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f ",i,a[i],b[i],c[i], as,bs,cs, factor);
#endif
		
        a[i]=as*factor-a[i];
        b[i]=bs*factor-b[i];
        c[i]=cs*factor-c[i];
		
#if defined(DEBUG)
        printf("out: %7.5f %7.5f\n",a[i],b[i],c[i]);
#endif
		
    }
}


//-------------------------------------------------------------------------------------------------------
// OverrelaxRedBulk
//-------------------------------------------------------------------------------------------------------
__global__ void OverrelaxRedBulk(REAL *a, REAL *b, REAL *c,
                                 REAL *an, REAL *bn, REAL *cn,
                                 REAL *Jpx, REAL *Jpy, REAL *Jpz,
                                 REAL *Jmx, REAL *Jmy, REAL *Jmz) {
	
    unsigned int iy, iz;
    int diff;
    REAL as, bs, cs;
	
    int spx, smx, spy, smy, spz, smz;
	
    REAL factor;
	
    int i;
	
    const unsigned int tid = threadIdx.x;
	
    const unsigned int num_threads = gridDim.x*blockDim.x;
    int istart, iend;
	
    istart=blockDim.x*blockIdx.x+d_Lsdb2; /*start on second plane (by using d_Lsdb2)
                                            since first plane is for the boundary */
    iend=d_nspin-d_Lsdb2;
	
    for(i=istart+tid; i<iend; i+=num_threads) {
		
        iz=i>>d_sdbLsm1;
        iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
        //do the same thing if both $iz and $iy are even or odd
        diff=(iz&1)^(iy&1);
		
        spx=(i >> (-diff)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (diff-1));
        smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (-diff)) + (i>>(diff-1));
        spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
        smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
        spz=(i+d_Lsdb2) /* & d_maskL3 */;
        smz=(i-d_Lsdb2) /* & d_maskL3 */;
		
        as=an[spx]*Jpx[i]+
            an[spy]*Jpy[i]+
            an[spz]*Jpz[i]+
            an[smx]*Jmx[i]+
            an[smy]*Jmy[i]+
            an[smz]*Jmz[i];
		
        bs=bn[spx]*Jpx[i]+
            bn[spy]*Jpy[i]+
            bn[spz]*Jpz[i]+
            bn[smx]*Jmx[i]+
            bn[smy]*Jmy[i]+
            bn[smz]*Jmz[i];
		
        cs=cn[spx]*Jpx[i]+
            cn[spy]*Jpy[i]+
            cn[spz]*Jpz[i]+
            cn[smx]*Jmx[i]+
            cn[smy]*Jmy[i]+
            cn[smz]*Jmz[i];
		
        cs+=d_h;
		
        //Overrelaxation
		
        factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/(as*as+bs*bs+cs*cs);
		
#if defined(DEBUG)
        printf("Reds [%d]: %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f ",i,a[i],b[i],c[i], as,bs,cs, factor);
#endif
		
        a[i]=as*factor-a[i];
        b[i]=bs*factor-b[i];
        c[i]=cs*factor-c[i];
		
#if defined(DEBUG)
        printf("out: %7.5f %7.5f\n",a[i],b[i],c[i]);
#endif
		
		
    }
	
}


//-------------------------------------------------------------------------------------------------------
// OverrelaxBlueBnd
//-------------------------------------------------------------------------------------------------------
__global__ void OverrelaxBlueBnd(REAL *a, REAL *b, REAL *c,
                                 REAL *an, REAL *bn, REAL *cn,
                                 REAL *Jpx, REAL *Jpy, REAL *Jpz,
                                 REAL *Jmx, REAL *Jmy, REAL *Jmz) {
	
    unsigned int iy, iz;
    int diff;
    REAL as, bs, cs;
	
    int spx, smx, spy, smy, spz, smz;
	
    REAL factor;
	
    int i;
	
    const unsigned int tid = threadIdx.x;
	
    const unsigned int num_threads = gridDim.x*blockDim.x;
    int istart;
	
    istart=blockDim.x*blockIdx.x; /* does first plane (by using d_Lsdb2) */
	
    for(i=istart+tid; i<d_Lsdb2; i+=num_threads) {
		
        iz=i>>d_sdbLsm1;
        iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
        //do the same thing if both $iz and $iy are even or odd
        diff=(iz&1)^(iy&1);
		
        spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
        smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
        spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
        smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
        spz=(i+d_Lsdb2) /* & d_maskL3 */;
        smz=(i-d_Lsdb2) /* & d_maskL3 */;
		
        as=an[spx]*Jpx[spx]+
            an[spy]*Jpy[spy]+
            an[spz]*Jpz[spz]+
            an[smx]*Jmx[smx]+
            an[smy]*Jmy[smy]+
            an[smz]*Jmz[smz];
		
        bs=bn[spx]*Jpx[spx]+
            bn[spy]*Jpy[spy]+
            bn[spz]*Jpz[spz]+
            bn[smx]*Jmx[smx]+
            bn[smy]*Jmy[smy]+
            bn[smz]*Jmz[smz];
		
        cs=cn[spx]*Jpx[spx]+
            cn[spy]*Jpy[spy]+
            cn[spz]*Jpz[spz]+
            cn[smx]*Jmx[smx]+
            cn[smy]*Jmy[smy]+
            cn[smz]*Jmz[smz];
		
        cs+=d_h;
		
        //Overrelaxation
		
        factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/
            (as*as+bs*bs+cs*cs);
#if defined(DEBUG)
        printf("Blues [%d]: %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f ",i,a[i],b[i],c[i], as,bs,cs, factor);
#endif
        a[i]=as*factor-a[i];
        b[i]=bs*factor-b[i];
        c[i]=cs*factor-c[i];
#if defined(DEBUG)
        printf("out: %7.5f %7.5f\n",a[i],b[i],c[i]);
#endif
		
    }
	
    istart=(d_nspin-d_Lsdb2)+blockDim.x*blockIdx.x;
	
    for(i=istart+tid; i<d_nspin; i+=num_threads) {
        iz=i>>d_sdbLsm1;
        iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
        //do the same thing if both $iz and $iy are even or odd
        diff=(iz&1)^(iy&1);
		
        spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
        smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
        spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
        smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
        spz=(i+d_Lsdb2) /* & d_maskL3 */;
        smz=(i-d_Lsdb2) /* & d_maskL3 */;
		
        as=an[spx]*Jpx[spx]+
            an[spy]*Jpy[spy]+
            an[spz]*Jpz[spz]+
            an[smx]*Jmx[smx]+
            an[smy]*Jmy[smy]+
            an[smz]*Jmz[smz];
		
        bs=bn[spx]*Jpx[spx]+
            bn[spy]*Jpy[spy]+
            bn[spz]*Jpz[spz]+
            bn[smx]*Jmx[smx]+
            bn[smy]*Jmy[smy]+
            bn[smz]*Jmz[smz];
		
        cs=cn[spx]*Jpx[spx]+
            cn[spy]*Jpy[spy]+
            cn[spz]*Jpz[spz]+
            cn[smx]*Jmx[smx]+
            cn[smy]*Jmy[smy]+
            cn[smz]*Jmz[smz];
		
        cs+=d_h;
		
        //Overrelaxation
		
        factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/
            (as*as+bs*bs+cs*cs);
#if defined(DEBUG)
        printf("Blues [%d]: %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f ",i,a[i],b[i],c[i], as,bs,cs, factor);
#endif
        a[i]=as*factor-a[i];
        b[i]=bs*factor-b[i];
        c[i]=cs*factor-c[i];
#if defined(DEBUG)
        printf("out: %7.5f %7.5f\n",a[i],b[i],c[i]);
#endif
    }
	
}


//-------------------------------------------------------------------------------------------------------
// OverrelaxBlueBulk
//-------------------------------------------------------------------------------------------------------
__global__ void OverrelaxBlueBulk(REAL *a, REAL *b, REAL *c,
                                  REAL *an, REAL *bn, REAL *cn,
                                  REAL *Jpx, REAL *Jpy, REAL *Jpz,
                                  REAL *Jmx, REAL *Jmy, REAL *Jmz) {
	
    unsigned int iy, iz;
    int diff;
    REAL as, bs, cs;
	
    int spx, smx, spy, smy, spz, smz;
	
    REAL factor;
	
    int i;
	
    const unsigned int tid = threadIdx.x;
	
    const unsigned int num_threads = gridDim.x*blockDim.x;
    int istart, iend;
	
    istart=blockDim.x*blockIdx.x; /*start on second plane (by using d_Lsdb2)
                                    since first plane is for the boundary */
    iend=d_nspin-d_Lsdb2;
	
    for(i=istart+tid; i<iend; i+=num_threads) {
		
        iz=i>>d_sdbLsm1;
        iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
        //do the same thing if both $iz and $iy are even or odd
        diff=(iz&1)^(iy&1);
		
        spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
        smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
        spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
        smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
        spz=(i+d_Lsdb2) /* & d_maskL3 */;
        smz=(i-d_Lsdb2) /* & d_maskL3 */;
		
        as=an[spx]*Jpx[spx]+
            an[spy]*Jpy[spy]+
            an[spz]*Jpz[spz]+
            an[smx]*Jmx[smx]+
            an[smy]*Jmy[smy]+
            an[smz]*Jmz[smz];
		
        bs=bn[spx]*Jpx[spx]+
            bn[spy]*Jpy[spy]+
            bn[spz]*Jpz[spz]+
            bn[smx]*Jmx[smx]+
            bn[smy]*Jmy[smy]+
            bn[smz]*Jmz[smz];
		
        cs=cn[spx]*Jpx[spx]+
            cn[spy]*Jpy[spy]+
            cn[spz]*Jpz[spz]+
            cn[smx]*Jmx[smx]+
            cn[smy]*Jmy[smy]+
            cn[smz]*Jmz[smz];
		
        cs+=d_h;
		
        //Overrelaxation
		
        factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/
            (as*as+bs*bs+cs*cs);
#if defined(DEBUG)
        printf("Blues [%d]: %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f ",i,a[i],b[i],c[i], as,bs,cs, factor);
#endif
        a[i]=as*factor-a[i];
        b[i]=bs*factor-b[i];
        c[i]=cs*factor-c[i];
#if defined(DEBUG)
        printf("out: %7.5f %7.5f\n",a[i],b[i],c[i]);
#endif
		
    }
	
}


#endif


//-------------------------------------------------------------------------------------------------------
// ComputeEnergy
//-------------------------------------------------------------------------------------------------------
__global__ void ComputeEnergy(REAL *a, REAL *b, REAL *c,
                              REAL *an, REAL *bn, REAL *cn,
                              REAL *Jpx, REAL *Jpy, REAL *Jpz,
                              REAL *Jmx, REAL *Jmy, REAL *Jmz,
                              REAL *energy,
                              int even) {
	
    unsigned int iy, iz;
    int diff;
    int spx, smx, spy, smy, spz, smz;
	
    int i;
	
    const unsigned int tid = threadIdx.x;
	
    const unsigned int num_threads = gridDim.x*blockDim.x;
    int istart;
	
    istart=blockDim.x*blockIdx.x; /*start on second plane (by using d_Lsdb2)
                                    since first plane is for the boundary */
	
	
    for(i=istart+tid; i<d_nspin; i+=num_threads) {
		
        iz=i>>d_sdbLsm1;
        iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
        //do the same thing if both $iz and $iy are even or odd
        diff=(((iz&1)^(iy&1))^even);
		
        spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
        smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
        spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
        smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
        spz=(i+d_Lsdb2) /* & d_maskL3 */;
        smz=(i-d_Lsdb2) /* & d_maskL3 */;
		
		
        energy[i]=a[i]*(an[spx]*Jpx[i]+
                        an[spy]*Jpy[i]+
                        an[spz]*Jpz[i]+
                        an[smx]*Jmx[i]+
                        an[smy]*Jmy[i]+
						an[smz]*Jmz[i])+
            b[i]*(bn[spx]*Jpx[i]+
                  bn[spy]*Jpy[i]+
                  bn[spz]*Jpz[i]+
                  bn[smx]*Jmx[i]+
                  bn[smy]*Jmy[i]+
                  bn[smz]*Jmz[i])+
            c[i]*(cn[spx]*Jpx[i]+
                  cn[spy]*Jpy[i]+
                  cn[spz]*Jpz[i]+
                  cn[smx]*Jmx[i]+
                  cn[smy]*Jmy[i]+
                  cn[smz]*Jmz[i]);
#if 0
		cuPrintf("i=%d, spx=%d, spy=%d, spz=%d, smx=%d, smy=%d, smz=%d, ", i, spx, spy, spz, smx, smy, smz);
		cuPrintf("a=%f, b=%f, c=%f, anspx=%f, anspy=%f, anspz=%f, ansmx=%f, ansmy=%f, ansmz=%f, ", a[i], b[i], c[i], an[spx],an[spy],an[spz], an[smx],an[smy],an[smz]);
		
		cuPrintf("bnspx=%f, bnspy=%f, bnspz=%f, bnsmx=%f, bnsmy=%f, bnsmz=%f, ", bn[spx],bn[spy],bn[spz], bn[smx],bn[smy],bn[smz]);
		cuPrintf("cnspx=%f, cnspy=%f, cnspz=%f, cnsmx=%f, cnsmy=%f, cnsmz=%f, ", cn[spx],cn[spy],cn[spz], cn[smx],cn[smy],cn[smz]);
		cuPrintf("Jpx=%f, Jpy=%f, Jpz=%f, Jmx=%f, Jmy=%f, Jmz=%f\n" , Jpx[i],Jpy[i],Jpz[i], Jmx[i],Jmy[smy],Jmz[i]);
#endif
	}
	
}


//-------------------------------------------------------------------------------------------------------
// ComputeHEnergy
//-------------------------------------------------------------------------------------------------------
void ComputeHEnergy(REAL *a, REAL *b, REAL *c,
                    REAL *an, REAL *bn, REAL *cn,
                    REAL *Jpx, REAL *Jpy, REAL *Jpz,
                    REAL *Jmx, REAL *Jmy, REAL *Jmz,
                    REALENE *energy,
                    int even) {
	
	unsigned int iy, iz;
	int diff;
	int spx, smx, spy, smy, spz, smz;
	REAL e;
	
	int i;
	
	for(i=0; i<h_nspin; i++) {
		iz=i>>h_sdbLsm1;
		iy=(i - (iz<<h_sdbLsm1))>>h_sdbLm1;
		
		//do the same thing if both $iz and $iy are even or odd
		diff=(((iz&1)^(iy&1))^even);
		
		spx=(i >> (diff-1)) + (((i-(i&h_maskL)) + ((i+1)&h_maskL)) >> (-diff));
		smx=(((i-(i&h_maskL)) + ((i-1)&h_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&h_maskL2))+((i+h_Ldb2)&h_maskL2);
		smy=(i-(i&h_maskL2))+((i-h_Ldb2)&h_maskL2);
		spz=(i+h_Lsdb2) & h_maskL3;
		smz=(i-h_Lsdb2) & h_maskL3;
		
		e=a[i]*(an[spx]*Jpx[i]+
				an[spy]*Jpy[i]+
				an[spz]*Jpz[i]+
				an[smx]*Jmx[i]+
				an[smy]*Jmy[i]+
				an[smz]*Jmz[i])+
            b[i]*(bn[spx]*Jpx[i]+
                  bn[spy]*Jpy[i]+
                  bn[spz]*Jpz[i]+
                  bn[smx]*Jmx[i]+
                  bn[smy]*Jmy[i]+
                  bn[smz]*Jmz[i])+
            c[i]*(cn[spx]*Jpx[i]+
                  cn[spy]*Jpy[i]+
                  cn[spz]*Jpz[i]+
                  cn[smx]*Jmx[i]+
                  cn[smy]*Jmy[i]+
                  cn[smz]*Jmz[i]);
		energy[0]+=e;
		
	}
	
}


//-------------------------------------------------------------------------------------------------------
// ComputeRandomGPU
//-------------------------------------------------------------------------------------------------------
void ComputeRandomGPU(int numBlocksRand, int numThreadsPerBlockRand, int mode, REAL* array){
	
	seedRand(0);
	
	MY_CUDA_CHECK( cudaThreadSynchronize() );
	
	dim3 dimGrid(numBlocksRand,1,1);
	
	dim3 dimBlock(numThreadsPerBlockRand,1,1);
	
	RandomGPU<<<dimGrid, numThreadsPerBlockRand>>>(array,mode);
	
	//cutilCheckMsg("RandomGPU() execution failed\n");
	
	MY_CUDA_CHECK( cudaThreadSynchronize() );
	
}


//-------------------------------------------------------------------------------------------------------
// FreeRand
//-------------------------------------------------------------------------------------------------------
void FreeRand(void){
	
	MY_CUDA_CHECK( cudaFree(d_aRand) );
	MY_CUDA_CHECK( cudaFree(d_bRand) );
	MY_CUDA_CHECK( cudaFree(d_cRand) );
	MY_CUDA_CHECK( cudaFree(d_Rand) );
	
}


//-------------------------------------------------------------------------------------------------------
// initRand
//-------------------------------------------------------------------------------------------------------
int initRand(int randomNumber){
	
	
	int n_per_rng;
	// int randNum; DA DISCOMMENTARE DOPO PROVA MESSO GLOBAL
	
	n_per_rng = iAlignUp(iDivUp(randomNumber, MT_RNG_COUNT), 2);
	
	randNum = MT_RNG_COUNT * n_per_rng;
	
    //printf("randomNumber: %d n_per_rng: %d MT_RNG_COUNT: %d randNum: %d\n",randomNumber,n_per_rng,MT_RNG_COUNT,randNum);
  
	MY_CUDA_CHECK( cudaMalloc((void **)&d_aRand, randNum * sizeof(REAL)) );
	MY_CUDA_CHECK( cudaMalloc((void **)&d_bRand, randNum * sizeof(REAL)) );
	MY_CUDA_CHECK( cudaMalloc((void **)&d_cRand, randNum * sizeof(REAL)) );
	MY_CUDA_CHECK( cudaMalloc((void **)&d_Rand, randNum * sizeof(REAL)) );
	
	MY_CUDA_CHECK( cudaMemcpyToSymbol("NPerRng",&n_per_rng,sizeof(int),0,cudaMemcpyHostToDevice) );
	
	return randNum;
	
}


//-------------------------------------------------------------------------------------------------------
// seedRand
//-------------------------------------------------------------------------------------------------------
void seedRand(int seed){
	
	if(seed==0) seed=(int)(drand48()*1023); // UNA SCHIFEZZA
	
	static mt_struct_stripped h_MT[MT_RNG_COUNT];
	static mt_struct MT[MT_RNG_COUNT];
	
	FILE *fd = Fopen("./data/MersenneTwister.raw", "rb");
	
	for (int i = 0; i < MT_RNG_COUNT; i++){
		//Inline structure size for compatibility,
		//since pointer types are 8-byte on 64-bit systems (unused *state variable)
		Fread(MT + i, 16 /* sizeof(mt_struct) */ * sizeof(int), 1, fd);
		
	}
	
	fclose(fd);
	
	fd = Fopen("./data/MersenneTwister.dat", "rb");
	
	Fread(h_MT, sizeof(h_MT), 1, fd);
	
	fclose(fd);
	
	//Need to be thread-safe
	mt_struct_stripped *MT_temp = (mt_struct_stripped *)Malloc(MT_RNG_COUNT * sizeof(mt_struct_stripped));
	
	for(int i = 0; i < MT_RNG_COUNT; i++){
		MT_temp[i] = h_MT[i];
		MT_temp[i].seed = seed;
	}
	
	MY_CUDA_CHECK( cudaMemcpyToSymbol("ds_MT",MT_temp,sizeof(h_MT),0,cudaMemcpyHostToDevice) );
	
	free(MT_temp);
	
}


//-------------------------------------------------------------------------------------------------------
// printRand
//-------------------------------------------------------------------------------------------------------
void printRand(int randomNumber,REAL* array,char* namefile){
	
	REAL * h_Rand  = (REAL *)malloc(randomNumber * sizeof(REAL));
	
	MY_CUDA_CHECK( cudaMemcpy(h_Rand, array, randomNumber * sizeof(REAL), cudaMemcpyDeviceToHost) );
	
	FILE* out = Fopen(namefile, "w");
	
	for(int x=0;x<randomNumber;x++) { fprintf(out,"%f\n",h_Rand[x]); }
	
	fclose(out);
	
}


////////////////////////////////////////////////////////////////////////////////
// Write MT_RNG_COUNT vertical lanes of NPerRng random numbers to *d_Random.
// For coalesced global writes MT_RNG_COUNT should be a multiple of warp size.
// Initial states for each generator are the same, since the states are
// initialized from the global seed. In order to improve distribution properties
// on small NPerRng supply dedicated (local) seed to each twister.
// The local seeds, in their turn, can be extracted from global seed
// by means of any simple random number generator, like LCG.
////////////////////////////////////////////////////////////////////////////////
__global__ void RandomGPU( REAL *d_Random, int mode ){
	
	const int tid = blockDim.x * blockIdx.x + threadIdx.x;
	const int THREAD_N = blockDim.x * gridDim.x;
	
	int iState, iState1, iStateM, iOut;
	unsigned int mti, mti1, mtiM, x;
	unsigned int mt[MT_NN];
	
	REAL temp_rand;
	
	for(int iRng = tid; iRng < MT_RNG_COUNT; iRng += THREAD_N){
		//Load bit-vector Mersenne Twister parameters
		mt_struct_stripped config = ds_MT[iRng];
		
		//Initialize current state
		mt[0] = config.seed;
		for(iState = 1; iState < MT_NN; iState++)
            mt[iState] = (1812433253U * (mt[iState - 1] ^ (mt[iState - 1] >> 30)) + iState) & MT_WMASK;
		
		iState = 0;
		mti1 = mt[0];
		
		for(iOut = 0; iOut < NPerRng; iOut++){
			//iState1 = (iState +     1) % MT_NN
			//iStateM = (iState + MT_MM) % MT_NN
			iState1 = iState + 1;
			iStateM = iState + MT_MM;
			if(iState1 >= MT_NN) iState1 -= MT_NN;
			if(iStateM >= MT_NN) iStateM -= MT_NN;
			mti  = mti1;
			mti1 = mt[iState1];
			mtiM = mt[iStateM];
			
			x = (mti & MT_UMASK) | (mti1 & MT_LMASK);
			x = mtiM ^ (x >> 1) ^ ((x & 1) ? config.matrix_a : 0);
			mt[iState] = x;
			iState = iState1;
			
			//Tempering transformation
			x ^= (x >> MT_SHIFT0);
			x ^= (x << MT_SHIFTB) & config.mask_b;
			x ^= (x << MT_SHIFTC) & config.mask_c;
			x ^= (x >> MT_SHIFT1);
			
			temp_rand = ((REAL)x + 1.0f) / 4294967296.0f; //Convert to (0, 1] float
			
			if(mode==ZEROONE) { // write to global memory
				d_Random[iRng + iOut * MT_RNG_COUNT] = temp_rand;
			} else if(mode==ONEONE) {
				d_Random[iRng + iOut * MT_RNG_COUNT] = 2.0f*temp_rand-1.0f;
			}else {
				d_Random[iRng + iOut * MT_RNG_COUNT] = temp_rand - 0.5f;
			}
			
		}
	}
}


//-------------------------------------------------------------------------------------------------------
// HeatBath
//-------------------------------------------------------------------------------------------------------
#if defined(SINGLEHBKERNEL)
__global__ void HeatBath(REAL *a, REAL *b, REAL *c,
                         REAL *an, REAL *bn, REAL *cn,
                         REAL *Jpx, REAL *Jpy, REAL *Jpz,
                         REAL *Jmx, REAL *Jmy, REAL *Jmz,
                         REAL *atry, REAL *btry, REAL *ctry,
                         REAL *random, int even) {
	
	unsigned int iy, iz;
	
	int diff;
	REAL as, bs, cs;
	REAL anew, bnew, cnew;
	REAL hk, bhk, projection, pl, invnorm;
	
	REAL invexpbhk2;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart;
	
	istart=blockDim.x*blockIdx.x;
	
	
	for(i=istart+tid; i<d_nspin; i+=num_threads) {
		
		anew=atry[i];
		bnew=btry[i];
		pl=anew*anew+bnew*bnew;
		if(pl>d_const_zerotwentyfive) {
			continue;
		}
		cnew=ctry[i];
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		//do the same thing if both $iz and $iy are even or odd
		diff=(((iz&1)^(iy&1))^even);
		
		spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2) & d_maskL3;
		smz=(i-d_Lsdb2) & d_maskL3;
		
		as=an[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
            an[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
            an[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
            an[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
            an[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
            an[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];
		
		bs=bn[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
            bn[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
            bn[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
            bn[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
            bn[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
            bn[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];
		
		cs=cn[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
            cn[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
            cn[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
            cn[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
            cn[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
            cn[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];
		
		cs+=d_h;
		
		hk=SQRT(as*as+bs*bs+cs*cs);
		bhk=d_Beta*hk;
		
		invexpbhk2=EXP(-d_const_two*bhk);
		projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
		if(projection<-d_const_one) { projection=-d_const_one; }
		
		invnorm=SQRT((d_const_one-cnew*cnew)/pl);
		
		anew*=invnorm;
		bnew*=invnorm;
		
		factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);
		
		anew = anew-factor*as;
		bnew = bnew-factor*bs;
		cnew = cnew-factor*cs;
		
		factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
		projection/=hk;
		a[i]=as*projection+anew*factor;
		b[i]=bs*projection+bnew*factor;
		c[i]=cs*projection+cnew*factor;
		
	}
	
}

#else


//-------------------------------------------------------------------------------------------------------
// HeatBathRedBnd
//-------------------------------------------------------------------------------------------------------
__global__ void HeatBathRedBnd(REAL *a, REAL *b, REAL *c,
							   REAL *an, REAL *bn, REAL *cn,
							   REAL *Jpx, REAL *Jpy, REAL *Jpz,
							   REAL *Jmx, REAL *Jmy, REAL *Jmz,
							   REAL *atry, REAL *btry, REAL *ctry,
							   REAL *random) {
	
	unsigned int iy, iz;
	
	int diff;
	REAL as, bs, cs;
	REAL anew, bnew, cnew;
	REAL hk, bhk, projection, pl, invnorm;
	
	REAL invexpbhk2;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart;
	
	istart=blockDim.x*blockIdx.x;
	
	
	for(i=istart+tid; i<d_Lsdb2; i+=num_threads) {
		anew=atry[i];
		bnew=btry[i];
		pl=anew*anew+bnew*bnew;
		if(pl>d_const_zerotwentyfive) {
			continue;
		}
		cnew=ctry[i];
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (-diff)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (diff-1));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (-diff)) + (i>>(diff-1));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2);// & d_maskL3;
		smz=(i-d_Lsdb2);// & d_maskL3;
		
		as=an[spx]*Jpx[i]+
            an[spy]*Jpy[i]+
            an[spz]*Jpz[i]+
            an[smx]*Jmx[i]+
            an[smy]*Jmy[i]+
            an[smz]*Jmz[i];
		
		bs=bn[spx]*Jpx[i]+
            bn[spy]*Jpy[i]+
            bn[spz]*Jpz[i]+
            bn[smx]*Jmx[i]+
            bn[smy]*Jmy[i]+
            bn[smz]*Jmz[i];
		
		cs=cn[spx]*Jpx[i]+
            cn[spy]*Jpy[i]+
            cn[spz]*Jpz[i]+
            cn[smx]*Jmx[i]+
            cn[smy]*Jmy[i]+
            cn[smz]*Jmz[i];
		
		cs+=d_h;
		
		hk=SQRT(as*as+bs*bs+cs*cs);
		bhk=d_Beta*hk;
		
		invexpbhk2=EXP(-d_const_two*bhk);
		projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
		if(projection<-d_const_one) { projection=-d_const_one; }
		
		invnorm=SQRT((d_const_one-cnew*cnew)/pl);
		
		anew*=invnorm;
		bnew*=invnorm;
		
		factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);
		
		anew = anew-factor*as;
		bnew = bnew-factor*bs;
		cnew = cnew-factor*cs;
		
		factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
		projection/=hk;
		a[i]=as*projection+anew*factor;
		b[i]=bs*projection+bnew*factor;
		c[i]=cs*projection+cnew*factor;
		
	}
	
	istart=(d_nspin-d_Lsdb2)+blockDim.x*blockIdx.x;
	
	
    for(i=istart+tid; i<d_nspin; i+=num_threads) {
		
		anew=atry[i];
		bnew=btry[i];
		pl=anew*anew+bnew*bnew;
		if(pl>d_const_zerotwentyfive) {
			continue;
		}
		cnew=ctry[i];
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (-diff)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (diff-1));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (-diff)) + (i>>(diff-1));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2);// & d_maskL3;
		smz=(i-d_Lsdb2);// & d_maskL3;
		
		as=an[spx]*Jpx[i]+
            an[spy]*Jpy[i]+
            an[spz]*Jpz[i]+
            an[smx]*Jmx[i]+
            an[smy]*Jmy[i]+
            an[smz]*Jmz[i];
		
		bs=bn[spx]*Jpx[i]+
            bn[spy]*Jpy[i]+
            bn[spz]*Jpz[i]+
            bn[smx]*Jmx[i]+
            bn[smy]*Jmy[i]+
            bn[smz]*Jmz[i];
		
		cs=cn[spx]*Jpx[i]+
            cn[spy]*Jpy[i]+
            cn[spz]*Jpz[i]+
            cn[smx]*Jmx[i]+
            cn[smy]*Jmy[i]+
            cn[smz]*Jmz[i];
		
		cs+=d_h;
		
		hk=SQRT(as*as+bs*bs+cs*cs);
		bhk=d_Beta*hk;
		
		invexpbhk2=EXP(-d_const_two*bhk);
		projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
		if(projection<-d_const_one) { projection=-d_const_one; }
		
		invnorm=SQRT((d_const_one-cnew*cnew)/pl);
		
		anew*=invnorm;
		bnew*=invnorm;
		
		factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);
		
		anew = anew-factor*as;
		bnew = bnew-factor*bs;
		cnew = cnew-factor*cs;
		
		factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
		projection/=hk;
		a[i]=as*projection+anew*factor;
		b[i]=bs*projection+bnew*factor;
		c[i]=cs*projection+cnew*factor;
		
	}
	
}


//-------------------------------------------------------------------------------------------------------
// HeatBathRedBulk
//-------------------------------------------------------------------------------------------------------
__global__ void HeatBathRedBulk(REAL *a, REAL *b, REAL *c,
								REAL *an, REAL *bn, REAL *cn,
								REAL *Jpx, REAL *Jpy, REAL *Jpz,
								REAL *Jmx, REAL *Jmy, REAL *Jmz,
								REAL *atry, REAL *btry, REAL *ctry,
								REAL *random) {
	
	unsigned int iy, iz;
	
	int diff;
	REAL as, bs, cs;
	REAL anew, bnew, cnew;
	REAL hk, bhk, projection, pl, invnorm;
	
	REAL invexpbhk2;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart, iend;
	
	istart=blockDim.x*blockIdx.x+d_Lsdb2; //start on second plane (by using d_Lsdb2)
										  // since first plane is for the boundary 
	iend=d_nspin-d_Lsdb2;
	
	for(i=istart+tid; i<iend; i+=num_threads) {
    
		anew=atry[i];
		bnew=btry[i];
		pl=anew*anew+bnew*bnew;
		if(pl>d_const_zerotwentyfive) {
			continue;
		}
		cnew=ctry[i];

		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (-diff)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (diff-1));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (-diff)) + (i>>(diff-1));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2);// & d_maskL3;
		smz=(i-d_Lsdb2);// & d_maskL3;
		
		as=an[spx]*Jpx[i]+
            an[spy]*Jpy[i]+
            an[spz]*Jpz[i]+
            an[smx]*Jmx[i]+
            an[smy]*Jmy[i]+
            an[smz]*Jmz[i];
		
		bs=bn[spx]*Jpx[i]+
            bn[spy]*Jpy[i]+
            bn[spz]*Jpz[i]+
            bn[smx]*Jmx[i]+
            bn[smy]*Jmy[i]+
            bn[smz]*Jmz[i];
		
		cs=cn[spx]*Jpx[i]+
            cn[spy]*Jpy[i]+
            cn[spz]*Jpz[i]+
            cn[smx]*Jmx[i]+
            cn[smy]*Jmy[i]+
            cn[smz]*Jmz[i];
		
		cs+=d_h;
		
		hk=SQRT(as*as+bs*bs+cs*cs);
		bhk=d_Beta*hk;
		
		invexpbhk2=EXP(-d_const_two*bhk);

		projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;

		if(projection<-d_const_one) { projection=-d_const_one; }
		
		invnorm=SQRT((d_const_one-cnew*cnew)/pl);
		
		anew*=invnorm;
		bnew*=invnorm;
		
		factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);
		
		anew = anew-factor*as;
		bnew = bnew-factor*bs;
		cnew = cnew-factor*cs;
		
		factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
		projection/=hk;
		a[i]=as*projection+anew*factor;
		b[i]=bs*projection+bnew*factor;
		c[i]=cs*projection+cnew*factor;
		
		
	}

}



//-------------------------------------------------------------------------------------------------------
// HeatBathBlueBnd
//-------------------------------------------------------------------------------------------------------
__global__ void HeatBathBlueBnd(REAL *a, REAL *b, REAL *c,
								REAL *an, REAL *bn, REAL *cn,
								REAL *Jpx, REAL *Jpy, REAL *Jpz,
								REAL *Jmx, REAL *Jmy, REAL *Jmz,
								REAL *atry, REAL *btry, REAL *ctry,
								REAL *random) {
	
	unsigned int iy, iz;
	
	int diff;
	REAL as, bs, cs;
	REAL anew, bnew, cnew;
	REAL hk, bhk, projection, pl, invnorm;
	
	REAL invexpbhk2;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart;
	
	istart=blockDim.x*blockIdx.x;
	
	
	for(i=istart+tid; i<d_Lsdb2; i+=num_threads) {
		
		anew=atry[i];
		bnew=btry[i];
		pl=anew*anew+bnew*bnew;
		if(pl>d_const_zerotwentyfive) {
			continue;
		}
		cnew=ctry[i];
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2);// & d_maskL3;
		smz=(i-d_Lsdb2);// & d_maskL3;
		
		as=an[spx]*Jpx[spx]+
            an[spy]*Jpy[spy]+
            an[spz]*Jpz[spz]+
            an[smx]*Jmx[smx]+
            an[smy]*Jmy[smy]+
            an[smz]*Jmz[smz];
		
		bs=bn[spx]*Jpx[spx]+
            bn[spy]*Jpy[spy]+
            bn[spz]*Jpz[spz]+
            bn[smx]*Jmx[smx]+
            bn[smy]*Jmy[smy]+
            bn[smz]*Jmz[smz];
		
		
		cs=cn[spx]*Jpx[spx]+
            cn[spy]*Jpy[spy]+
            cn[spz]*Jpz[spz]+
            cn[smx]*Jmx[smx]+
            cn[smy]*Jmy[smy]+
            cn[smz]*Jmz[smz];
		
		cs+=d_h;
		
		hk=SQRT(as*as+bs*bs+cs*cs);
		bhk=d_Beta*hk;
		
		invexpbhk2=EXP(-d_const_two*bhk);
		projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
		if(projection<-d_const_one) { projection=-d_const_one; }
		
		invnorm=SQRT((d_const_one-cnew*cnew)/pl);
		
		anew*=invnorm;
		bnew*=invnorm;
		
		factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);
		
		anew = anew-factor*as;
		bnew = bnew-factor*bs;
		cnew = cnew-factor*cs;
		
		factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
		projection/=hk;
		a[i]=as*projection+anew*factor;
		b[i]=bs*projection+bnew*factor;
		c[i]=cs*projection+cnew*factor;
		
	}
	
	istart=(d_nspin-d_Lsdb2)+blockDim.x*blockIdx.x;
	
	
	for(i=istart+tid; i<d_nspin; i+=num_threads) {
		
		anew=atry[i];
		bnew=btry[i];
		pl=anew*anew+bnew*bnew;
		if(pl>d_const_zerotwentyfive) {
			continue;
		}
		cnew=ctry[i];
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2);// & d_maskL3;
		smz=(i-d_Lsdb2);// & d_maskL3;
		
		as=an[spx]*Jpx[spx]+
            an[spy]*Jpy[spy]+
            an[spz]*Jpz[spz]+
            an[smx]*Jmx[smx]+
            an[smy]*Jmy[smy]+
            an[smz]*Jmz[smz];
		
		bs=bn[spx]*Jpx[spx]+
            bn[spy]*Jpy[spy]+
            bn[spz]*Jpz[spz]+
            bn[smx]*Jmx[smx]+
            bn[smy]*Jmy[smy]+
            bn[smz]*Jmz[smz];
		
		
		cs=cn[spx]*Jpx[spx]+
            cn[spy]*Jpy[spy]+
            cn[spz]*Jpz[spz]+
            cn[smx]*Jmx[smx]+
            cn[smy]*Jmy[smy]+
            cn[smz]*Jmz[smz];
		
		cs+=d_h;
		
		hk=SQRT(as*as+bs*bs+cs*cs);
		bhk=d_Beta*hk;
		
		invexpbhk2=EXP(-d_const_two*bhk);
		projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
		if(projection<-d_const_one) { projection=-d_const_one; }
		
		invnorm=SQRT((d_const_one-cnew*cnew)/pl);
		
		anew*=invnorm;
		bnew*=invnorm;
		
		factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);
		
		anew = anew-factor*as;
		bnew = bnew-factor*bs;
		cnew = cnew-factor*cs;
		
		factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
		projection/=hk;
		a[i]=as*projection+anew*factor;
		b[i]=bs*projection+bnew*factor;
		c[i]=cs*projection+cnew*factor;
		
	}
	
	
}



//-------------------------------------------------------------------------------------------------------
// HeatBathBlueBulk
//-------------------------------------------------------------------------------------------------------
__global__ void HeatBathBlueBulk(REAL *a, REAL *b, REAL *c,
								 REAL *an, REAL *bn, REAL *cn,
								 REAL *Jpx, REAL *Jpy, REAL *Jpz,
								 REAL *Jmx, REAL *Jmy, REAL *Jmz,
								 REAL *atry, REAL *btry, REAL *ctry,
								 REAL *random) {
	
	
	unsigned int iy, iz;
	
	int diff;
	REAL as, bs, cs;
	REAL anew, bnew, cnew;
	REAL hk, bhk, projection, pl, invnorm;
	
	REAL invexpbhk2;
	
	int spx, smx, spy, smy, spz, smz;
	
	REAL factor;
	
	int i;
	
	const unsigned int tid = threadIdx.x;
	
	const unsigned int num_threads = gridDim.x*blockDim.x;
	int istart, iend;
	
	istart=blockDim.x*blockIdx.x;
	//istart=blockDim.x*blockIdx.x+d_Lsdb2; //start on second plane (by using d_Lsdb2)
    // since first plane is for the boundary 

	iend=d_nspin-d_Lsdb2;
	
	for(i=istart+tid; i<iend; i+=num_threads) {
		
		anew=atry[i];
		bnew=btry[i];
		pl=anew*anew+bnew*bnew;
		if(pl>d_const_zerotwentyfive) {
			continue;
		}
		cnew=ctry[i];
		iz=i>>d_sdbLsm1;
		iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;
		
		diff=(iz&1)^(iy&1);
		
		spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
		smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
		spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
		smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
		spz=(i+d_Lsdb2);// & d_maskL3;
		smz=(i-d_Lsdb2);// & d_maskL3;
		
		as=an[spx]*Jpx[spx]+
            an[spy]*Jpy[spy]+
            an[spz]*Jpz[spz]+
            an[smx]*Jmx[smx]+
            an[smy]*Jmy[smy]+
            an[smz]*Jmz[smz];
		
		bs=bn[spx]*Jpx[spx]+
            bn[spy]*Jpy[spy]+
            bn[spz]*Jpz[spz]+
            bn[smx]*Jmx[smx]+
            bn[smy]*Jmy[smy]+
            bn[smz]*Jmz[smz];
		
		
		cs=cn[spx]*Jpx[spx]+
            cn[spy]*Jpy[spy]+
            cn[spz]*Jpz[spz]+
            cn[smx]*Jmx[smx]+
            cn[smy]*Jmy[smy]+
            cn[smz]*Jmz[smz];
		
		cs+=d_h;
		
		hk=SQRT(as*as+bs*bs+cs*cs);
		bhk=d_Beta*hk;
		
		invexpbhk2=EXP(-d_const_two*bhk);
		projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
		if(projection<-d_const_one) { projection=-d_const_one; }
		
		invnorm=SQRT((d_const_one-cnew*cnew)/pl);
		
		anew*=invnorm;
		bnew*=invnorm;
		
		factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);
		
		anew = anew-factor*as;
		bnew = bnew-factor*bs;
		cnew = cnew-factor*cs;
		
		factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
		projection/=hk;
		/*a[i]=as*projection+anew*factor;
          b[i]=bs*projection+bnew*factor;
          c[i]=cs*projection+cnew*factor;*/
		
	}
	
}


#endif

struct mpiexch *SetUpMpi(int id, int np, int side) {
	struct mpiexch *p;
	p=(struct mpiexch*)Malloc(sizeof(struct mpiexch));
	MY_CUDA_CHECK( cudaMallocHost((void **)&p->upsend,side*side*sizeof(MYREAL)*NCOMP) );
	MY_CUDA_CHECK( cudaMallocHost((void **)&p->uprecv,side*side*sizeof(MYREAL)*NCOMP) );
	MY_CUDA_CHECK( cudaMallocHost((void **)&p->downsend,side*side*sizeof(MYREAL)*NCOMP) );
	MY_CUDA_CHECK( cudaMallocHost((void **)&p->downrecv,side*side*sizeof(MYREAL)*NCOMP) );
	p->up=id+1<np?id+1:0;
	p->down=id-1<0?np-1:id-1;
	p->nbyte=side*side*sizeof(MYREAL)*NCOMP;
	p->post=FALSE;
	return p;
}

/*
 * Local variables:
 *  mode: c++
 *  c-indent-level: 4
 *  c-basic-offset: 4
 *  tab-width: 4
 *  indent-tabs-mode: nil
 * End:
 */
